Summary

We obtain the IRFs for some of the best models (with respect to AIC/BIC/Shapiro-Wilk/Jarque-Bera/Ljung-Box). We use either estimated values for the static shock transmission matrix \(B\) or impose the long-run restriction suggested by Blanchard and Quah (1989).

In addition, we generate the IRFs of Blanchard and Quah (1989) and GMR.

Preliminaries

Load and plot data

tt = readRDS(paste0(params$PATH, params$JOBID, "/",
                    "e_residualcheck_end.rds")) %>% arrange(rk_bic)

DATASET = readRDS(params$PATH_DATASET)
dygraphs::dygraph(DATASET)
DATASET = DATASET %>% as.matrix()
DIM_OUT = dim(DATASET)[2]
N_OBS = dim(DATASET)[1]

Define some functions for IRFs

# Permute static shock transmission matrix B and change signs

permute_chgsign = function(irf_array, 
                           perm = rep(1, dim(irf_array)[1]), 
                           sign = rep(1, dim(irf_array)[1])){
  
  dim_out = dim(irf_array)[1]
  
  perm_mat = diag(dim_out)
  perm_mat = perm_mat[,perm]
  
  sign = diag(sign)
  
  ll = map(1:dim(irf_array)[3], ~ irf_array[,,.x] %*% perm_mat %*% sign)
  
  irf_array = array(0, dim = c(dim_out, dim_out, dim(irf_array)[3]))
  for (ix in 1:dim(irf_array)[3]){
    irf_array[,,ix] = ll[[ix]] 
  }
  
  return(irf_array)
}

# Plot IRF

plot_irf = function(irf_array){
  plot(pseries(irf_array %>% polm(), lag.max = 40))
}


# Calculate the cumulative sum of the IRF coefficients for a given element 

irf_cumsum = function(irf_array, el_mat){
  
  stopifnot("*el_mat* should be a matrix with two columns where each row contains the indices for which the cumulative sum should be calculated" = 
            dim(el_mat)[2] == 2)
  
  for (ix in (1:nrow(el_mat))){
    irf_array[el_mat[ix,1], el_mat[ix,2], ] = irf_array[el_mat[ix,1], el_mat[ix,2], ] %>% cumsum()
  }
  
  return(irf_array)
}


# Rotate the static shock transmission matrix such that there is no long-run effect

rotate_longrun = function(irf_array){

  dim_out = dim(irf_array)[1]
  
  stopifnot("Number of outputs must be equal to 2" = dim_out == 2)
  
  k1 = freqresp(armamod(lmfd(polm(diag(dim_out)), polm(irf_array)), sigma_L = diag(dim_out)), n.f = 1)$frr
  lq = lq_decomposition(matrix(c(Re(k1)), nrow = dim_out))
  Q = lq$q

  jj = map(1:dim(irf_array)[3], ~ irf_array[,,.x] %*% t(Q))
  
  irf_array = array(0, dim = c(dim_out, dim_out,dim(irf_array)[3]))
  for (ix in 1:dim(irf_array)[3]){
    irf_array[,,ix] = jj[[ix]] 
  }

return(irf_array)
}

Model statistics

tt %>%
  select(-nr, -p_plus_q, -rk_mle, -contains("value"), -contains("cov"), -contains("pval")) %>% 
  filter(normality_flag == 0, lb_flag == 0) %>% 
  dtable()
## This version of bslib is designed to work with rmarkdown version 2.7 or higher.

Blanchard Quah AR(8)

Function that calculates the cumulative IRF for some variables of choice. Here, we are interested in log(GNP) and not its change.

Estimate a VAR(8) model with the svars package.

vars_bq = VAR(DATASET, p = 8, type = "none")
vars_irf_bq = irf(vars_bq, n.ahead = 100)

Transform it to rationalmatrices/RLDM style and create a first plot.

irf_bq = array(0, dim = c(2,2,101))
irf_bq[,1,] = t(vars_irf_bq$irf$rGDPgrowth_demeaned)
irf_bq[,2,] = t(vars_irf_bq$irf$unemp_detrended)

plot(pseries(polm(irf_bq), lag.max = 100))

Calculate the long-run response, LQ decompose it, rotate the static shock impact matrix, transform the transfer function s.t. the demand shock is transitory.

k1 = matrix(c(colSums(vars_irf_bq$irf$rGDPgrowth_demeaned),
              colSums(vars_irf_bq$irf$unemp_detrended)),
            nrow = 2)
lq = lq_decomposition(k1)
Q = lq$q
 
ii = map(1:101, ~ irf_bq[,,.x] %*% t(Q))

irf_bq_lr = array(0, dim = c(2,2,101))
for (ix in 1:101){
  irf_bq_lr[,,ix] = ii[[ix]] 
}
rm(k1, Q, ii, irf_bq)

Check the result : transfer function evaluated at \(1\) should be triangular.

apply(irf_bq_lr, c(1,2), sum) %>% zapsmall()
##           [,1]     [,2]
## [1,]  0.622691 0.000000
## [2,] -0.161077 5.558043

Plot the rotated IRF

plot(pseries(irf_bq_lr %>% polm(), lag.max = 40))

irf_bq_lr_cumsum = irf_cumsum(irf_bq_lr, matrix(c(1,1,1,2), nrow = 2))

Plot the rotated IRF where the response to (log) GNP is cumulated.

plot(pseries(irf_bq_lr_cumsum %>% polm(), lag.max = 40))

Permute columns and change sign such that the plots align with the ones in Blanchard and Quah.

jj = map(1:101, ~ irf_bq_lr_cumsum[,,.x] %*% matrix(c(0,-1,1,0), nrow = 2))

irf_bq_lr_cumsum_permuted = array(0, dim = c(2,2,101))
for (ix in 1:101){
  irf_bq_lr_cumsum_permuted[,,ix] = jj[[ix]] 
}

plot(pseries(irf_bq_lr_cumsum_permuted %>% polm(), lag.max = 40))

GMR VARMA(4,1, k = 1)

if(!file.exists("../local_data/gmr_results_bq.rds")){
  load("~/r_projects/gmr_ssvarma/results/BQ/results.MLE.4lags.res")
  write_rds(res.estim.MLE, "../local_data/gmr_results_bq.rds")
  gmr_results_bq = res.estim.MLE
} else{
  gmr_results_bq = read_rds("../local_data/gmr_results_bq.rds")
}

Convert GMR results to rationalmatrices/RLDM style.

ar_gmr = gmr_results_bq$Phi
ma_gmr = gmr_results_bq$Theta
B_gmr = gmr_results_bq$C

polm_ar = array(c(c(diag(2)),
                  -c(ar_gmr)),
                dim = c(2,2,5)) %>% polm()

polm_ma = array(c(c(diag(2)),
                  -c(ma_gmr)),
                dim = c(2,2,2)) %>% polm()

B = gmr_results_bq$C

lmfd_obj = lmfd(polm_ar, polm_ma %r% B)

Check poles and zeros

polm_ar %>% zeroes() %>% abs()
## [1] 1.192445 1.192445 1.324240 1.447863 1.447863 2.216027 2.216027 4.102746
polm_ma %>% zeroes() %>% abs()
## [1] 0.565828 2.519183
irf_gmr = pseries(lmfd_obj, lag.max = 40) %>% unclass()
plot_irf(irf_gmr)

irf_gmr_cumsum = irf_cumsum(irf_gmr, 
                           el_mat = matrix(c(1,1,1,2), nrow = 2))

plot_irf(irf_gmr_cumsum)

Permute columns such that results align

irf_gmr_cumsum_permuted = permute_chgsign(irf_gmr_cumsum, perm = c(2,1))
plot_irf(irf_gmr_cumsum_permuted)

BF VARMA(1,2, kappa = 1, k = 0)

Extract appropriate inter-valued parameters.

tt %>% 
  filter(p==1 & q==2 & kappa==1 & k==0) %>% 
  select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_12_2

Create an ARMA-WHF model:

bf_params_12_2 = tt_bf_12_2$params_deep_final %>% .[[1]]
bf_tmpl_12_2 = tt_bf_12_2$tmpl %>% .[[1]]

write_rds(bf_params_12_2, "../local_data/p_whf/best_mod_rob/bf_params_12_2.rds")
write_rds(bf_tmpl_12_2, "../local_data/p_whf/best_mod_rob/bf_tmpl_12_2.rds")


bf_armawhfmod_12_2 = fill_tmpl_whf_rev(bf_params_12_2, bf_tmpl_12_2)

Create IRF, cumulate GNP, and plot:

bf_irf_12_2 = irf_whf(bf_params_12_2, bf_tmpl_12_2, n_lags = 40) %>% unclass()
bf_irf_12_2_cumsum = bf_irf_12_2 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_12_2_cumsum %>% plot_irf()

bf_irf_12_2_cumsum_permuted = bf_irf_12_2_cumsum

Rotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.

bf_irf_12_2_lr = bf_irf_12_2 %>% rotate_longrun()
bf_irf_12_2_lr_cumsum = bf_irf_12_2_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_12_2_lr_cumsum %>% plot_irf()

bf_irf_12_2_lr_cumsum_permuted = bf_irf_12_2_lr_cumsum %>% permute_chgsign(perm = c(2,1),
                                                                           sign = c(-1,1))
bf_irf_12_2_lr_cumsum_permuted %>% plot_irf()

BF VARMA(1,3, kappa = 1, k = 0)

Extract appropriate inter-valued parameters.

tt %>% 
  filter(p==1 & q==3 & kappa==1 & k==0) %>% 
  select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_13_2

Create an ARMA-WHF model:

bf_params_13_2 = tt_bf_13_2$params_deep_final %>% .[[1]]
bf_tmpl_13_2 = tt_bf_13_2$tmpl %>% .[[1]]

bf_armawhfmod_13_2 = fill_tmpl_whf_rev(bf_params_13_2, bf_tmpl_13_2)

Create IRF, cumulate GNP, and plot:

bf_irf_13_2 = irf_whf(bf_params_13_2, bf_tmpl_13_2, n_lags = 40) %>% unclass()
bf_irf_13_2_cumsum = bf_irf_13_2 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_13_2_cumsum %>% plot_irf()

bf_irf_13_2_cumsum_permuted = bf_irf_13_2_cumsum

Rotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.

bf_irf_13_2_lr = bf_irf_13_2 %>% rotate_longrun()
bf_irf_13_2_lr_cumsum = bf_irf_13_2_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_13_2_lr_cumsum %>% plot_irf()

bf_irf_13_2_lr_cumsum_permuted = bf_irf_13_2_lr_cumsum %>% permute_chgsign(perm = c(2,1),
                                                                           sign = c(-1,1))
bf_irf_13_2_lr_cumsum_permuted %>% plot_irf()

BF VARMA(3,1, kappa = 1, k = 0)

Extract appropriate inter-valued parameters.

tt %>% 
  filter(p==3 & q==1 & kappa==1 & k==0) %>% 
  select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_31_2

Create an ARMA-WHF model:

bf_params_31_2 = tt_bf_31_2$params_deep_final %>% .[[1]]
bf_tmpl_31_2 = tt_bf_31_2$tmpl %>% .[[1]]

bf_armawhfmod_31_2 = fill_tmpl_whf_rev(bf_params_31_2, bf_tmpl_31_2)

Create IRF, cumulate GNP, and plot:

bf_irf_31_2 = irf_whf(bf_params_31_2, bf_tmpl_31_2, n_lags = 40) %>% unclass()
bf_irf_31_2_cumsum = bf_irf_31_2 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_31_2_cumsum %>% plot_irf()

Sign-permute such that shocks are labelled consistently.

bf_irf_31_2_cumsum_permuted = bf_irf_31_2_cumsum %>% permute_chgsign(perm = c(2,1),
                                                                     sign = c(1,1))
bf_irf_31_2_cumsum_permuted %>% plot_irf()

Rotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.

bf_irf_31_2_lr = bf_irf_31_2 %>% rotate_longrun()
bf_irf_31_2_lr_cumsum = bf_irf_31_2_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_31_2_lr_cumsum %>% plot_irf()

bf_irf_31_2_lr_cumsum_permuted = bf_irf_31_2_lr_cumsum %>% permute_chgsign(perm = c(2,1),
                                                                           sign = c(-1,1))
bf_irf_31_2_lr_cumsum_permuted %>% plot_irf()

BF VARMA(1,1, kappa = 1)

Extract appropriate inter-valued parameters.

tt %>% 
  filter(p==1 & q==1 & kappa==1 & k==0) %>% 
  select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_11_2

Create an ARMA-WHF model:

bf_params_11_2 = tt_bf_11_2$params_deep_final %>% .[[1]]
bf_tmpl_11_2 = tt_bf_11_2$tmpl %>% .[[1]]

bf_armawhfmod_11_2 = fill_tmpl_whf_rev(bf_params_11_2, bf_tmpl_11_2)

Create IRF, cumulate GNP, and plot:

bf_irf_11_2 = irf_whf(bf_params_11_2, bf_tmpl_11_2, n_lags = 40) %>% unclass()
bf_irf_11_2_cumsum = bf_irf_11_2 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_11_2_cumsum %>% plot_irf()

Sign-permute such that shocks are labelled consistently.

bf_irf_11_2_cumsum_permuted = bf_irf_11_2_cumsum %>% permute_chgsign(perm = c(2,1),
                                                                     sign = c(-1,-1))
bf_irf_11_2_cumsum_permuted %>% plot_irf()

Rotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.

bf_irf_11_2_lr = bf_irf_11_2 %>% rotate_longrun()
bf_irf_11_2_lr_cumsum = bf_irf_11_2_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_11_2_lr_cumsum %>% plot_irf()

bf_irf_11_2_lr_cumsum_permuted = bf_irf_11_2_lr_cumsum %>% permute_chgsign(perm = c(2,1),
                                                                           sign = c(-1,1))
bf_irf_11_2_lr_cumsum_permuted %>% plot_irf()

BF VARMA(4,1, kappa = 0, k = 1): GMR’s model

Extract appropriate inter-valued parameters.

tt %>% 
  filter(p==4 & q==1 & kappa==0 & k==1) %>% 
  select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_41_1

Create an ARMA-WHF model:

bf_params_41_1 = tt_bf_41_1$params_deep_final %>% .[[1]]
bf_tmpl_41_1 = tt_bf_41_1$tmpl %>% .[[1]]

bf_armawhfmod_41_1 = fill_tmpl_whf_rev(bf_params_41_1, bf_tmpl_41_1)

Create IRF, cumulate GNP, and plot:

bf_irf_41_1 = irf_whf(bf_params_41_1, bf_tmpl_41_1, n_lags = 40) %>% unclass()
bf_irf_41_1_cumsum = bf_irf_41_1 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_41_1_cumsum %>% plot_irf()

Sign-permute such that shocks are labelled consistently.

bf_irf_41_1_cumsum_permuted = bf_irf_41_1_cumsum %>% permute_chgsign(perm = c(2,1))
bf_irf_41_1_cumsum_permuted %>% plot_irf()

Rotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.

bf_irf_41_1_lr = bf_irf_41_1 %>% rotate_longrun()
bf_irf_41_1_lr_cumsum = bf_irf_41_1_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_41_1_lr_cumsum %>% plot_irf()

bf_irf_41_1_lr_cumsum_permuted = bf_irf_41_1_lr_cumsum %>% permute_chgsign(perm = c(2,1),
                                                                           sign = c(-1,1))
bf_irf_41_1_lr_cumsum_permuted %>% plot_irf()

Summary

list_mods = tibble(BQ_AR8 = irf_bq_lr_cumsum_permuted[,,1:41] %>% list(),
                   GMR_ARMA41_1 = irf_gmr_cumsum_permuted %>% list(),
                   BF_ARMA11_2 = bf_irf_11_2_cumsum_permuted %>% list(),
                   BF_ARMA11_2_lr = bf_irf_11_2_lr_cumsum_permuted %>% list(),
                   BF_ARMA12_2 = bf_irf_12_2_cumsum_permuted %>% list(),
                   BF_ARMA12_2_lr = bf_irf_12_2_lr_cumsum_permuted %>% list(),
                   BF_ARMA13_2 = bf_irf_13_2_cumsum_permuted %>% list(),
                   BF_ARMA13_2_lr = bf_irf_13_2_lr_cumsum_permuted %>% list(),
                   BF_ARMA31_2 = bf_irf_31_2_cumsum_permuted %>% list(),
                   BF_ARMA31_2_lr = bf_irf_31_2_lr_cumsum_permuted %>% list(),
                   BF_ARMA41_1 = bf_irf_41_1_cumsum_permuted %>% list(),
                   BF_ARMA41_1_lr = bf_irf_41_1_lr_cumsum_permuted %>% list())

list_mods %>% pivot_longer(everything()) %>% 
  mutate(Demand2GNP = map(value, ~.x[1,1,]),
         Demand2Unempl = map(value, ~.x[2,1,]),
         Supply2GNP = map(value, ~.x[1,2,]),
         Supply2Unempl = map(value, ~.x[2,2,])) %>% 
  select(-value) %>% 
  pivot_longer(c("Demand2GNP", "Demand2Unempl", "Supply2GNP", "Supply2Unempl"), 
               names_to = "Response_Type", values_to = "Impact") %>% 
  separate(Response_Type, into = c("Shock", "Output"), sep = "2") %>% 
  unnest_longer(Impact, indices_include = TRUE) %>% 
  rename(Lag = Impact_id,
         Model = name) %>% 
  mutate(Lag = Lag-1) -> tibble_irf
tibble_irf %>% 
  # filter(Model != "BF_ARMA11_2") %>% 
  ggplot(aes(x = Lag, y = Impact, color = Model)) +
  geom_line() + geom_point() + facet_grid(Output~Shock) -> ply

plotly::ggplotly(ply)
knitr::knit_exit()